% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function [Xi,Ti,Ui,cell_type] = el_resample(grid,base,U,varargin)
% [Xi,Ti,Ui] = el_resample(grid,base,U)
% [Xi,Ti,Ui] = el_resample(grid,base,U,k)
% [Xi,Ti,Ui] = el_resample(grid,base,U,k,l)
%
% Resample the solution U (stored as two dimensional array) on a local
% mesh for each element. The two optional arguments k and l are the
% global and the local element degrees. If both are absent, base.k is
% used; if l is absent, l=k is used.
%
% Note: the local mesh is defined in "ResampleMesh.octxt"

 if(nargin<4)
   k = base.k;
 else
   k = varargin{1};
 end
 if(nargin<5)
   l = k;
 else
   l = varargin{2};
 end
   
 locmeshes = load("ResampleMesh.octxt","meshes");
 if(k>size(locmeshes.meshes,2))
   error(['Global degree ',num2str(k),' not supported.']);
 endif
 if(l>size(locmeshes.meshes,3))
   error(['Local degree ',num2str(l),' not supported.']);
 endif

 locmesh = locmeshes.meshes(grid.d,k,l);
 p = locmesh.p;
 t = locmesh.t;

 % 1) evaluate the basis functions on the local mesh
 np = size(p,2);
 nt = size(t,2);
 PHI = zeros(base.pk,np);
 for i=1:base.pk
   PHI(i,:) = ev_pol(base.p_s{i},p);
 end

 % 2) resample all the elements
 Xi = zeros(grid.d,np,grid.ne);
 Ti = t; % the connectivity is always the same
 Ui = zeros(       np,grid.ne);
 idx = ones(1,np);
 for ie=1:grid.ne
   Xi(:,:,ie) = grid.e(ie).b * p + grid.e(ie).x0(:,idx);
   Ui(:,ie)   = PHI' * U(:,ie);
 end

 % 3) Paraview specific information
 if(nargout>3)
   cell_type = locmesh.PV_cell_type;
 endif

return

